\documentclass{article}
\usepackage{lmodern}
\usepackage[a4paper,scale=0.8]{geometry}
\usepackage{listings}
\usepackage[dvipsnames]{xcolor}
\lstset{
    language=C++,
    basicstyle=\ttfamily,
    breaklines=true,
    keywordstyle=\bfseries\color{NavyBlue}, 
    commentstyle=\itshape\color{black!50!white},
    stringstyle=\bfseries\color{PineGreen!90!black},
    columns=flexible,
    numbers=left,
    numbersep=2em,
    numberstyle=\footnotesize,
    frame=single,
    framesep=1em
}

\usepackage{pgf-umlcd}

\usepackage{ctex}
\usepackage{amsmath}
\usepackage{graphicx}
\usepackage{float}
\usepackage{datetime}

\title{数值分析作业1}
\author{陈乐瑶  3210103710}
\date{\today}

\begin{document}
\maketitle

\part{三种方法的类实现}
\begin{tikzpicture}
    \begin{class}[text width=5cm]{EquationSolver}{5,0}
        \attribute{+eps: double}
        \attribute{+max: double}
        \operation{+EquationSolver(e,m)}
        \operation{$\sim$EquationSolver()}
        \operation[0]{+solve(condition):double}
    \end{class}
    \begin{class}[text width=4cm]{Bisection}{0,-5}
        \inherit{EquationSolver}
        \attribute{+del: double}
        \operation{+Bisection(e,d,m)}
        \operation{$\sim$Bisection()}
        \operation{+solve(condition):double}
    \end{class}
    \begin{class}[text width=4cm]{Newton}{5,-5}
        \inherit{EquationSolver}
        \operation{+Newton(e,m)}
        \operation{$\sim$Newton()}
        \operation{+solve(condition):double}
    \end{class}
    \begin{class}[text width=4cm]{Secant}{10,-5}
        \inherit{EquationSolver}
        \attribute{+del: double}
        \operation{+Secant(e,d,m)}
        \operation{$\sim$Secant()}
        \operation{+solve(condition):double}
    \end{class}
    \begin{class}[text width=5cm]{Function}{10,-10}
        \operation[0]{+get\_value(x):double}
        \operation{+get\_derivation(x):double}
    \end{class}
    \begin{class}[text width=5cm]{InitInput}{1,-9}
        \attribute{+a,b: double}
        \attribute{+x0,x1: double}
        \attribute{+f: Function}
        \operation{+setBisectionInput(a,b):void}
        \operation{+setNewtonInput(x0):void}
        \operation{+setSecantInput(x0,x1):void}
    \end{class}
    \unidirectionalAssociation{Bisection}{}{}{InitInput}
    \unidirectionalAssociation{Newton}{}{}{InitInput}
    \unidirectionalAssociation{Secant}{}{}{InitInput}
    \unidirectionalAssociation{InitInput}{}{}{Function}
\end{tikzpicture}

\begin{itemize}
    \item Function类的继承类在不同实际函数测试cpp文件中给出
    \item 三种方法的solve函数基于讲义算法实现，在此不一一列出
    \item $f(x),f^{'}(x)$分别由Function类中的$f.get\_value(x)$和$f.get\_derivation(x)$实现
\end{itemize}


\part{三种方法的测试}
在三种方法的测试（B、C、D题）中，默认设置$\epsilon=10^{-8}$，$\delta=10^{-8}$，$M=30$。
分别用程序测试，得到近似根及该根的函数值，并以函数的几何作图为参考比照。
\section{B:bisection method}
以第一个函数$ \frac{1}{x}-\tan(x) $为例给出部分代码：
\begin{lstlisting}
class func1: public Function
{
    public:
        double get_value(double x){
            if (x == 0){
                cout << "Error: divided by zero!" << endl;
                return 0;
            }else if (cos(x) == 0){
                cout << "Error: 'tan(x)' does not exist!" <<endl;
                return 0;
            }else
                return 1.0/x-tan(x);
        }
}f1;
int main()
{
    cout << "Bisection Method Test: "<<endl;
    double e = 1e-8;    //set the epsilon
    double d = 1e-8;    //set the delta
    double m = 30;      //set the M
    
    //build a bisection solver with the above parameter
    Bisection bis_solver(e,d,m);

    //set the initial input
    InitInput input1(f1);
    input1.setBisectionInput(e, M_PI_2);

    //solve the equation
    double x1 = bis_solver.solve(input1);
    cout << "the approximate root of f1: x1 = "<< x1 <<endl;
    cout << "the value of the root: f1(x1) = "<< f1.get_value(x1) <<endl;
    cout <<endl;
}
\end{lstlisting}

需要注意的是，（1）、（2）在给出的区间端点存在无定义的情况，因此可对端点$\pm \epsilon$。\\
另一方面，（4）在给定区间无零点，且存在渐近线无值，因此迭代结果不是正确的零点。

运行结果及函数几何对照图如下：
\begin{figure}[htbp]
    \includegraphics[width=15cm]{./pic/B_result.jpg}
    \caption{Bisection test result}
\end{figure}
\begin{figure}[htbp]
    \begin{minipage}[t]{0.24\textwidth}
        \includegraphics[width=3cm]{./pic/B1.jpg}
        \caption{$ \frac{1}{x}-\tan\left(x\right) $}
    \end{minipage}
    \begin{minipage}[t]{0.24\textwidth}
        \includegraphics[width=3cm]{./pic/B2.jpg}
        \caption{$ \frac{1}{x}-2^{x} $}
    \end{minipage}
    \begin{minipage}[t]{0.24\textwidth}
        \includegraphics[width=3cm]{./pic/B3.jpg}
        \caption{$ 2^{-x}+e^{x}+2\cos(x)-6 $}
    \end{minipage}
    \begin{minipage}[t]{0.24\textwidth}
        \includegraphics[width=3cm]{./pic/B4.jpg}
        \caption{$ \frac{x^{3}+4x^{2}+3x+5}{2x^{3}-9x^{2}+18x-2} $}
    \end{minipage}
\end{figure}

\section{C:Newton's method}
\begin{lstlisting}
class func: public Function
{
    public:
        double get_value(double x){
            if (cos(x) == 0){
                cout << "Error: 'tan(x)' does not exist!" <<endl;
                return 0;
            }else
                return x-tan(x);
        }
        /* the actual derivation is replaced by limit value
        double get_derivation(double x){
            if (cos(x) == 0){
                cout << "Error: divided by zero!" <<endl;
                return 0;
            }else
                return 1.0-1.0/(cos(x)*cos(x));
        }*/
}f;
int main()
{
    cout << "Newton's Method Test: "<<endl;
    double e = 1e-8;    //set the epsilon
    double m = 30;      //set the M
    
    //build a Newton solver with the above parameter
    Newton new_solver(e, m);

    double x1 = 4.5, x2 = 7.7;  //set the two initial x
    InitInput input1(f), input2(f);//set the initial input
    input1.setNewtonInput(x1);
    double r1 = new_solver.solve(input1);//solve the equation
    cout << "the approximate root of f near 4.5: r1 = "<< r1 << endl;
    cout << "the value of the root: f(r1) = "<< f.get_value(r1) <<endl;
    cout << endl;

    return 0;
}
\end{lstlisting}

运行结果和几何对照图如下：
\begin{figure}[htbp]
    \includegraphics[width=15cm]{./pic/C_result.jpg}
    \caption{Newton test result}
\end{figure}
\begin{figure}[htbp]
    \begin{minipage}[t]{0.24\textwidth}
        \includegraphics[width=3cm]{./pic/c1.jpg}
        \caption{$ x-\tan(x) $}
    \end{minipage}
\end{figure}

\section{D:secant method}
以第一个函数$ \sin(\frac{x}{2})-1 $为例给出部分代码：
\begin{lstlisting}
class func1: public Function
{
    public:
        double get_value(double x){
            return sin(x/2.0)-1;
        }
}f1;
int main()
{
    cout << "Secant Method Test: "<<endl;
    double e = 1e-8;    //set the epsilon
    double d = 1e-8;    //set the delta
    double m = 30;      //set the M
    
    //build a secant solver with the above parameter
    Secant sec_solver(e,d,m);
 
    //set the initial input
    InitInput input1(f1);
    input1.setSecantInput(0, M_PI_2-e);

    //solve the equation
    double x1 = sec_solver.solve(input1);
    cout << "the approximate root of f1: x1 = "<< x1 <<endl;
    cout << "the value of the root: f1(x1) = "<< f1.get_value(x1) <<endl;
    cout << endl;

    return 0;
}
\end{lstlisting}

运行结果和几何对照图如下：
\begin{figure}[htbp]
    \includegraphics[width=15cm]{./pic/D_result.jpg}
    \caption{Secant test result}
\end{figure}
\begin{figure}[htbp]
    \begin{minipage}[t]{0.25\textwidth}
        \includegraphics[width=3cm]{./pic/D1.jpg}
        \caption{$ \sin(\frac{x}{2})-1 $}
    \end{minipage}
    \begin{minipage}[t]{0.25\textwidth}
        \includegraphics[width=3cm]{./pic/D2.jpg}
        \caption{$ e^{x}-\tan(x) $}
    \end{minipage}
    \begin{minipage}[t]{0.25\textwidth}
        \includegraphics[width=3cm]{./pic/D3.jpg}
        \caption{$ x^{3}-12x^{2}+3x+1 $}
    \end{minipage}
\end{figure}

通过更改x0,x1的初始值，可以得到不同的根值。
由函数的几何作图可以发现，该函数的零点有多个，因此不同的初始值可以迭代得到不同的根值。

\part{三种方法的问题应用}
\section{E:trough problem}
该问题可以抽象为：
\begin{equation}
    V = f(h) = L\bigl(0.5 \pi r^{2} - r^{2}\arcsin\frac{h}{r} - h{(r^{2} - h^{2})}^{\frac{1}{2}}\bigr) = 12.43ft^{3}
\end{equation}
其中，$ L=10ft, r=1ft$，误差需要控制在$ 0.01ft $。

\begin{lstlisting}
    class func: public Function
    {
        public:
            double L = 10.0;
            double r = 1.0;
            double V_standard = 12.4;  
            double get_value(double x){ 
                double V = L*(0.5*M_PI*r*r - r*r*asin(x/r) 
                            - x*sqrtf(r*r-x*x));
                return V-V_standard;
            }
            /* the actual derivation is replaced by limit value
            double get_derivation(double x){
                double V_der = -2*L*sqrtf(r*r-x*x);
                return V_der;
            }*/
    }f;
    int main()
    {
        double e = 0.0001;  //set the epsilon(smaller than the tolerance)
        double d = 0.0001;  //set the delta(smaller than the tolerance)
        double m = 50;      //set the M
    
        //build three solver for each method
        Bisection bis_solver(e, d, m);
        Newton new_solver(e, m);
        Secant sec_solver(e, d, m);
    
        InitInput input1(f), input2(f), input3(f);
        input1.setBisectionInput(0, 1);
        double h1 = bis_solver.solve(input1);
        cout << "Bisection method: h1 = " << h1 << endl;

        input2.setNewtonInput(0.1);
        double h2 = new_solver.solve(input2);
        cout << "Newton's method: h2 = " << h2 << endl;

        input3.setSecantInput(0, 1);
        double h3 = sec_solver.solve(input3);
        cout << "Secant method: h3 = " << h3 << endl;

        return 0;
    }
\end{lstlisting}

运行结果如下：
\begin{figure}[htbp]
    \includegraphics[width=15cm]{./pic/E_result.jpg}
    \caption{trough problem result}
\end{figure}

因此，槽的水深h应为0.166ft

\section{F:vehicle problem}
该问题可以抽象为：
\begin{gather}
    f(\alpha) = A\sin\alpha\cos\alpha + B\sin^{2}\alpha - C\cos\alpha -E\sin\alpha = 0\\
    A = l\sin\beta_{1}, B = l\cos\beta_{1},\\
    C = (h+0.5D)\sin\beta_{1} - 0.5D\tan\beta_{1},\\
    E = (h+0.5D)\cos\beta_{1} - 0.5D
\end{gather}
其中，$ l=89in, h=49in, D = 55in, \beta_{1}=11.5^{\circ}$。

\begin{lstlisting}
class func: public Function
{
    public:
        double l,h,D,beta;
        double A = l*sin(beta);
        double B = l*cos(beta);
        double C = (h+0.5*D)*sin(beta)-0.5*D*tan(beta);
        double E = (h+0.5*D)*cos(beta)-0.5*D;
    
        func(double l, double h, double D, double beta):
            l(l), h(h), D(D), beta(beta){};
        double get_value(double x){
            double val = A*sin(x)*cos(x) + B*sin(x)*sin(x) 
                        - C*cos(x) - E*sin(x);
            return val;
        }
        /* the actual derivation is replaced by limit value
        double get_derivation(double x){
            double der = A*(cos(x)*cos(x)-sin(x)*sin(x)) 
                        + 2*B*sin(x)*cos(x) + C*sin(x) - E*cos(x);
            return der;
        }*/
};
int main()
{
    double e = 1e-8;    //set the epsilon
    double d = 1e-8;    //set the delta
    double m = 50;      //set the M
    
    //build two solver for Newton method and secant method
    Newton new_solver(e, m);
    Secant sec_solver(e, d, m);
    
    //problem(a)
    func f1(89, 49, 55, 11.5);
    InitInput input1(f1);
    input1.setNewtonInput(33.0);
    double alpha1 = new_solver.solve(input1);
    cout << "(a) alpha = " << alpha1 << " (with val " << f1.get_value(alpha1) << ")" << endl;

    //problem(b): D=30in
    func f2(89, 49, 30, 11.5);
    InitInput input2(f2);
    input2.setNewtonInput(33.0);
    double alpha2 = new_solver.solve(input2);
    cout << "(b) alpha = " << alpha2 << " (with val " << f2.get_value(alpha2) << ")" << endl;

    //problem(c): set x1 away from 33
    InitInput input3(f1);
    input3.setSecantInput(33.0, 31.0);
    double alpha3 = sec_solver.solve(input3);
    cout << "(c) x1 = 31, alpha = " << alpha3 << " (with val " << f1.get_value(alpha3) << ")" << endl;
    input3.setSecantInput(33.0, 30.0);
    alpha3 = sec_solver.solve(input3);
    cout << "(c) x1 = 30, alpha = " << alpha3 << " (with val " << f1.get_value(alpha3) << ")" << endl;
    input3.setSecantInput(33.0, 16.0);
    alpha3 = sec_solver.solve(input3);
    cout << "(c) x1 = 16, alpha = " << alpha3 << " (with val " << f1.get_value(alpha3) << ")" << endl;
    input3.setSecantInput(33.0, 3.0);
    alpha3 = sec_solver.solve(input3);
    cout << "(c) x1 = 3, alpha = " << alpha3 << " (with val " << f1.get_value(alpha3) << ")" << endl;

    return 0;
}
\end{lstlisting}

运行结果如下：
\begin{figure}[htbp]
    \includegraphics[width=15cm]{./pic/F_result.jpg}
    \caption{vehicle problem result}
\end{figure}

对（c），经过多次测试发现，初始值x1设为（0，90）内不同值时会得到差异很大的不同根值。\\
猜测是因为该函数有多个零解，且函数曲线变化剧烈，不同点的导数差异很大，导致secant method得到的结果差异较大。\\
通过函数作图证实:
\begin{figure}[htbp]
    \includegraphics[width=8cm]{./pic/F_extra.jpg}
    \caption{vehicle function}
\end{figure}

\end{document}
